Dynamical critical exponent of the Jaynes-Cummings-Hubbard model 
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An array of high-Q electromagnetic resonators coupled to qubits gives rise to the Jaynes-Cummings-Hubbard 
model describing a superfluid to Mott insulator transition of lattice polaritons. From mean-field and strong 
coupling expansions, the critical properties of the model are expected to be identical to the scalar Bose-Hubbard 
model. A recent Monte Carlo study of the superfluid density on the square lattice suggested that this does not 
hold for the fixed-density transition through the Mott lobe tip. Instead, mean-field behavior with a dynamical 
critical exponent z — 2 was found. We perform large-scale quantum Monte Carlo simulations to investigate 
the critical behavior of the superfluid density and the compressibility. We find z — 1 at the tip of the insulating 
lobe. Hence the transition falls in the 3D XY universality class, analogous to the Bose-Hubbard model. 

PACS numbers: 71.36.-hc, 73.43.Nq, 78.20.Bh, 42.50.Ct 



Introduction A remarkable success of theoretical physics 
is the development of a theory of phase transitions which pro- 
vides a unified description of apparently distinct systems in 
terms of a few different universality classes characterized by 
universal critical exponents. The universality of a transition 
is completely determined by the dimension of the system, the 
order parameter and the range of interactions. An interesting 
example is the Mott insulator to superfluid (MI-SF) transition 
of lattice bosons, first studied in the framework of the Bose- 
Hubbard model (BHM) Q). When driven by density fluctua- 
tions, the MI-SF transition is known to be in the universality 
class of the dilute Bose gas with a dynamical critical exponent 
z = 2 (U 121. However, when the density of bosons is kept 
fixed during the transition its universality class changes and 
corresponds to the (2+l)D XY model with z = 1 |Q][fl|4). 

The experimental realization of the BHM with ultra-cold 
atoms in optical lattices nearly a decade ago has opened 
up a fast growing and versatile field of condensed matter 
physics. More recently, the realization of Bose-Einstein con- 
densation of weakly interacting polaritons, i.e., quasiparticles 
that form when light strongly interacts with matter, in high- 
er cavities Q has triggered an immense interest in realizing 
condensed matter systems with photonic systems. A recent 
theoretical focus has been on the MI-SF transition of polari- 
tons 17HIS1. The Jaynes-Cummings-Hubbard model (JCHM) 
has been introduced to describe such a quantum phase tran- 
sition of light in an array of coupled quantum electrodynam- 
ics (QED) cavities, each containing a single photonic mode 
interacting with a two-level system (qubit) [7-9|. The com- 
petition between strong photon-qubit coupling, giving rise to 
an effective photonic repulsion (localization), and the photon 
hopping between cavities (derealization) leads to a phase di- 
agram featuring Mott lobes reminiscent of those of ultracold 
atoms in optical lattices as described by the BHM. The real- 
ization of the JCHM has been proposed in various solid-state 
or quantum-optical systems, e.g., with nitrogen-vacancy (NV) 
centers in diamond GH9), quantum dot excitons |17|, super- 
conducting qubits [ 13 1 and trapped ions 1 18 1. Device integra- 
tion, high tunability and individual addressability of each cav- 
ity make wide parameter regimes easily accessible. Cavity or 



circuit QED arrays thus constitute one of the most promising 
architectures for quantum information processing and offer 
the possibility to study fundamental questions of interacting 
quantum systems. For a recent review of many-body physics 
in coupled-cavity arrays see Ref. IT9ll . 

The phase diagram and excitation spectra of the JCHM 
have been accurately determined by analytical Q [T3l4T5ll as 
well as numerical methods ll8T fl2"l [161 l20ll . However, an in- 
triguing open question concerns the universality class of the 
fixed-density transition of the JCHM. Quantum Monte Carlo 
(QMC) calculations of the superfluid density [12] suggest that 
the universality class at the tip of the Mott lobe of the JCHM 
is different (namely mean-field like) from the BHM. On the 
contrary, strong-coupling expansion lfl4l [T5l and an effective 
action approach 1 13 1 show that the same critical theory applies 
to both models. The discrepancy with the QMC predictions is 
very surprising, especially since the analytical arguments for 
arriving at a critical theory are similar to the BHM and are 
well understood. One would therefore expect the same scal- 
ing behavior, e.g., of the superfluid density, for the BHM and 
the JCHM. 

In this Rapid Communication we resolve the controversy 
between analytical and numerical findings and present results 
from extensive QMC simulations on the two-dimensional 
(2D) square lattice. We go beyond previous studies in two 
significant ways: (i) by using much larger system sizes, and 
(ii) by studying both the superfluid density and the compress- 
ibility. 

Model The JCHM is defined by the Hamiltonian 

H = ^h?-t^(a\a j +1Lc.)-i i N, (1) 

« (ij) 
with the JC Hamiltonian for site i ET1 

h\ c = Wpdja- + Wqo+o-r + g (a+a,. + crT/at) , 

and the total number of polaritons N = J2i( a l a i + <7 t (J 7) 
(combined number of photons and qubit excitations), which 
is a conserved quantity and can be controlled via the chemical 
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FIG. 1: Phase diagram of the JCHM with zero detuning on the 2D 
square lattice as obtained from QMC simulations 1 12]. Shown is the 
first Mott lobe with filling factor n = 1. The arrow indicates the 
fixed-density transition through the tip of the lobe. 



potential fi. Here, w p (w q ) is the energy of photons (qubits). 
The spin operators a ± describe intrasite transitions between 
the two qubit levels induced by emission or absorption of a 
photon with rate g (light-matter coupling). The second term in 
([TJ describes photon transfer between nearest neighbor sites i 
and j with hopping amplitude t. The polariton density (filling 
factor) is denoted as n = (n) = (N)/L 2 , where L 2 is the 
number of lattice sites. We use g as the unit of energy, and 
consider uj p = to q = 1 (resonance condition). 

Similar to the BHM, the ground-state phase diagram of the 
JCHM consists of a series of Mott-insulating lobes with inte- 
ger filling factor n. The extent of the lowest Mott lobe with 
n = 1 for the 2D square lattice considered here is known 
from QMC simulations |[T2l . see Fig. [TJ The generic MI-SF 
transition in the 2D BHM caused by density fluctuations is 
known to be of the mean-field type with a dynamical criti- 
cal exponent z — 2 JT] [2j. However, when increasing the 
hopping t along a path that leads through the tip of the Mott 
lobe while keeping the polariton density constant, the transi- 
tion is driven by phase fluctuations. This special case belongs 
to the universality class of the (2+l)D XY model with z = 1 
UJB). This means that particle and hole excitations become 
symmetric near the tip of the lobe and that the critical field 
theory describing this transition is relativistic with equal scal- 
ing of space and time directions. The location of the lobe tip 
has been determined as fjt/g = 0.185(5), t e /g = 0.05201(10) 
|[T2l . For the 2D JCHM, previous QMC simulations have con- 
firmed z — 2 for the generic, density-driven transition. How- 
ever, the same QMC simulations also suggest that the special 
z = 1 point at the tip of the lobe in the BHM is absent in the 
JCHM. On the other hand, analytical calculations of the exci- 
tation spectra lfl4l[T5l and general symmetry arguments based 
on the gauge invariance of the theory [ 1 3 1 predict a special 
point with z = 1 at the lobe tip in the JCHM as well. 

Method In order to resolve this controversy and to deter- 
mine the dynamical critical exponent z, we compute the finite- 
size scaling behavior of the superfluid density p s and the com- 



pressibility k. Scaling relations for both quantities are known 
from Ref. [TJ. The superfluid density is related to fluctuations 
of the winding number W in QMC simulations via 11221 

Ps ~ (3DL D - 2 ■ ( ' 

The finite-size scaling form of p s reads 

ft = L a - D -*ft[(t-t c )L 1 ^,)9/L*], (3) 

where v denotes the correlation length exponent. Fixing the 
ratio a = (3/L z , the quantity 

L D - a+ *p.[(t-tc)X 1/v ,a] (4) 

depends only on the distance from the critical point, (t — 
t c )L x l h ', so that at t — t c curves for different L as a function 
of t intersect in a single point. This allows to determine the 
critical value t c for the MI-SF transition. Plotting L D+Z ~ 2 p s 
as a function of (t — t c )L l ' v should lead to a scaling collapse 
of results for different L. 

A second observable of interest is the compressibility k = 
dn/dfi, which, by the fluctuation-dissipation theorem, can 
also be expressed in terms of number fluctuations, i.e., 

k = (3 ((n 2 ) - (h) 2 ) , (5) 

with the scaling form 

K = L z - D K[(t-t c )L 1/u ,a] (6) 

implying that L d ~ z k should be independent of L at t c . 

For the calculation of p s and n, the Hamiltonian ([TJl is sim- 
ulated using world lines in the stochastic series expansion 
(SSE) representation [ 20 , 23 1 . This is the same method as pre- 
viously used by Zhao et al. |fl~2). We have used the ALPS 1.3 
implementation [23 1 of the SSE with directed loop updates 
E4l|25). In the vicinity of the MI-SF boundary of the first 
Mott lobe, allowing a maximum of six photons respectively 
polaritons per site is sufficient to eliminate any noticeable er- 
ror. We considered Lx L square lattices. The inverse temper- 
ature (3 was chosen as j3g = 2L for z = 1 and /3g = L 2 /4 for 
z = 2. We also performed QMC simulations using the worm 
algorithm |26| in the path integral representation following 
the implementation of Ref. E71 . No cutoff in the polariton 
number is needed in this case. The results on up to 64 x 64 
lattices (not shown) are in full agreement with our SSE data. 

The conclusion of z — 2 scaling for the whole MI-SF phase 
boundary in previous work was based on results for the super- 
fluid density on lattice sizes up to 22 x 22 lfT2ll . Here we 
present data for both the superfluid density and the compress- 
ibility, using much larger system sizes up to 40 x 40. The 
existence of a z = 1 critical point should also be visible in 
numerical simulations in a finite range of parameters around 
the lobe tip [12|. Hence, a grand-canonical algorithm with a 
suitably chosen chemical potential (here /i/g = 0.185 |[T2l ) 
can be used, as confirmed by the results below. 
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FIG. 2: (Color online) Scaling of the superfiuid density p s across the 
fixed-density transition at fi/g = 0.185 1121 using L x L square 
lattices and j3g = 2L. The intersect of Lp s for different lattice sizes 
L in a single point in panel (a) is evidence for a dynamical critical 
exponent z = 1, and defines the critical point at t c /g = 0.05242(1). 
(b) Scaling collapse using tjg = 0.05242 and v = 0.6715 0. 



Results First, we consider the scaling behavior of p s and k 
assuming z = 1. In that case fluctuations in space and (imag- 
inary) time are isotropic and the winding number in space 
(leading to a nonzero superfiuid stiffness) scales the same 
way as the winding number in imaginary time (leading to a 
nonzero compressibility). Results for z = 1 with the aspect 
ratio fig/L = const. = 2 are shown in Figs.|2]and[3] 

Figure [2| a) shows the rescaled superfiuid density p s L as a 
function of the hopping strength t/g for system sizes ranging 
from 20 x 20 to 40 x 40. The intersect of the curves leads to the 
estimate of the critical hopping strength t c /g — 0.05241(1). 
The previous estimate was t c /g = 0.05201(10) [12|. Out- 
value of t c , together with the correlation length exponent v = 
0.6715 (as found numerically for the BHM [4]), leads to a 
clean scaling collapse in Fig. |2|b). We also observe that we 
need bigger system sizes for the JCHM as compared to the 
BHM in order to see a clear scaling collapse. 

Figure [3] shows a similar analysis for the compressibility n. 
The value of t c from the intersect in Fig.|3ja) coincides within 
errorbars with the value obtained from the superfiuid density. 



FIG. 3: (Color online) Scaling of the compressibility k for the same 
parameters as in Fig. [2] The intersect in (a) is consistent with z = 1 
and occurs at the same value of the critical hopping strength t c /g = 
0.05242(1) as in Fig.|2] (b) Scaling collapse using t c /g = 0.05242 
and v = 0.6715 0). 



Again, we find a very clean scaling collapse in Fig.|3jb). Our 
results are thus fully consistent with z = 1. 

In Ref. [12 1, an apparent scaling collapse for the superfiuid 
density was found assuming z = 2. In Fig. [4] we therefore 
present our results for the case z — 2 with f3g/L 2 — const. = 
1 /4. We point out that we do not find any contradiction be- 
tween our numerical data and those of Ref. [12] when con- 
sidering the same system sizes. The curves for the superfiuid 
density, shown in the Fig. [4]^ a X seem to cross in a single point 
and one may be tempted to think that z = 2 applies equally 
well 11121 . However, differences clearly show up for larger 
system sizes. FigureQb) shows finite-size scaling corrections 
to the critical value of the hopping strength that go as 1/L 2 . 
The extrapolation of those intersection points yields a criti- 
cal value for the hopping strength that is within error bars the 
same as the one found for z = 1 scaling. The estimate of t c for 
smaller L matches the result of Ref. [ 12 1. These observations 
can be understood from Eq. d2}: The winding number, which 
is an integer, is in 2D given by (W 2 ) ~ p s /3, and hence the 
leading term in the finite size scaling for the superfiuid den- 
sity cannot distinguish between z = 1 and z = 2. The small 
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FIG. 4: (Color online) Scaling of the superfluid density p s and com- 
pressibility k across the fixed-density transition at \ij g = 0.185 1121 
using jig — L 2 /4. In (a), the intersect of curves for consecutive sys- 
tem sizes [shown in panel (b)] approaches the i c obtained in Figs. [2] 
and[3] Panel (c) shows the compressibility for the same parameters, 
with no intersection in the vicinity of t c . Panel (d) illustrates the 
intersect obtained when plotting kL 2 . 



subleading corrections we see in Fig.^b) (while we were un- 
able to resolve any such drifts in the z — 1 scenario within 
the resolution of our numerical data) therefore provide further 
evidence that z = 1 is correct. Using the full range of system 
sizes available here, no scaling collapse is achieved assuming 
z = 2. 

An even stronger distinction between z = 1 and z = 2 
can be made using results for the compressibility. To this end, 
we plot in Fig. B[b) nL° = n as appropriate for z = 2 — 
cf. Eq. ([6j. The absence of a crossing point on the scale of 
the plot suggests that there are strong corrections to the scal- 
ing assumption, making z = 2 very unlikely. In contrast, k 
scales as expected for the generic transition of the JCHM (not 
shown). 

Finally, if z = 1 is correct, relativistic invariance allows 
us to interpret the imaginary time axis as a spatial axis, and 
vice versa. Consequently, kL 2 (with f3 ~ L 2 scaling) should 
behave identically to pL 2 . This is shown in Fig. P[d). We 
therefore conclude that z = 1 must be true, and we see no 
physical reason to investigate other scenarios such as theories 
with a fractional z. 

Summary We have shown on the basis of finite-size scal- 
ing of the superfluid density and the compressibility that the 
fixed-density Mott-insulator-superfluid transition in the 2D 
Jaynes-Cummings-Hubbard model falls into the 3D XY uni- 
versality class. The critical behavior is thus identical to the 
corresponding transition in the Bose-Hubbard model. 
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